slope <- map_dfr(rigal_trends, ~.x %>%
  select(siteid, linear_slope),
.id = "response"
)

slope_df <- slope %>%
  pivot_wider(names_from = "response", values_from = "linear_slope") %>%
  mutate(across(starts_with("log_"), ~(exp(.x) -1) * 100))

1 Check trends

p <- map2(rigal_trends, names(rigal_trends),
  ~.x %>%
    ggplot(aes(x = linear_slope)) +
    geom_histogram() +
    labs(title = .y)
)
plot_grid(plotlist = p)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

no NA in slopes since fix rigal_classification:

summary_slope <- slope %>%
  group_by(response) %>%
  summarise(summ = list(enframe(summary_distribution(linear_slope, na.rm = TRUE)))) %>%
  unnest(cols = summ) %>%
  pivot_wider(names_from = "name", values_from = "value")
summary_slope %>%
  mutate(response = get_var_replacement()[response]) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "bordered", "hover")) %>%
  scroll_box(width = "100%", height = "400px")

response	min	1st_quart	median	2nd_quart	max	mean	sd	n
Appearance	-0.0559441	-0.0018271	0.0009691	0.0065140	0.0600000	0.0027393	0.0098864	2678
Chao	-0.0975818	-0.0105013	-0.0015503	0.0006336	0.1177589	-0.0055673	0.0155877	2678
Chao evenness	-1.0762433	-0.0197209	0.0000000	0.0214427	0.8698342	0.0000118	0.0687473	2678
Chao species richness	-1.0804000	-0.0363727	0.0002926	0.0473422	0.9908720	0.0073944	0.1387268	2678
Chao shannon	-0.7782364	-0.0270535	0.0007389	0.0362457	0.7822804	0.0043759	0.0822483	2678
Chao simpson	-0.6747701	-0.0241523	0.0007189	0.0304959	0.5440252	0.0028965	0.0691378	2678
Disappearance	-0.0486529	-0.0005100	0.0000000	0.0042135	0.0557576	0.0017718	0.0071868	2678
Evenness	-0.0830806	-0.0060849	0.0000000	0.0069975	0.0817728	0.0006221	0.0144629	2678
Turnover abundance	-0.0781734	-0.0115266	-0.0037510	0.0014597	0.0672796	-0.0051970	0.0135836	2678
Horn	-0.0900632	-0.0125663	-0.0022471	0.0006332	0.0958701	-0.0064480	0.0169154	2678
Jaccard	-0.0931698	-0.0118699	-0.0045654	0.0011438	0.0548771	-0.0056807	0.0130380	2678
Jaccard dissimilarity	-0.0672796	-0.0014597	0.0037510	0.0115266	0.0781734	0.0051970	0.0135836	2678
Log species richness	-0.2307581	-0.0094470	0.0004979	0.0139823	0.2389473	0.0025976	0.0262344	2678
Log total abundance	-0.3751466	-0.0284703	0.0049807	0.0409316	0.4436236	0.0060785	0.0656325	2678
Nestedness	-0.0830556	-0.0041642	0.0002881	0.0074839	0.1156349	0.0021108	0.0137269	2678
Shannon	-0.2007711	-0.0088851	0.0004650	0.0121768	0.1574411	0.0016449	0.0221596	2678
Simpson	-0.0829149	-0.0045179	0.0001071	0.0058510	0.0800458	0.0007642	0.0111222	2678
Species richness	-1.6727273	-0.0328546	0.0020269	0.0536209	1.5727273	0.0148247	0.1422215	2678
Total turnover	-0.0559441	-0.0013253	0.0032257	0.0097964	0.0621302	0.0045111	0.0115567	2678
Total abundance	-755.8336538	-1.8770744	0.1129827	2.2456287	398.4499936	1.7786740	31.5332673	2678
turnover	-0.1223776	0.0000000	0.0000000	0.0077107	0.0989899	0.0030863	0.0134089	2678

na_siteid_slope <- unique(slope[is.na(slope$linear_slope), ]$siteid)

sd_num <- analysis_dataset %>%
  group_by(siteid) %>%
  summarise(across(where(is.numeric), ~(sd(.x))))

sd_num %>%
  filter(siteid %in% na_siteid_slope) %>%
  summary(.)
#>     siteid            species_nb  total_abundance log_total_abundance log_species_nb      year        latitude     longitude      main_bas   chao_richness  chao_shannon  chao_simpson
#>  Length:0           Min.   : NA   Min.   : NA     Min.   : NA         Min.   : NA    Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA  
#>  Class :character   1st Qu.: NA   1st Qu.: NA     1st Qu.: NA         1st Qu.: NA    1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA  
#>  Mode  :character   Median : NA   Median : NA     Median : NA         Median : NA    Median : NA   Median : NA   Median : NA   Median : NA   Median : NA   Median : NA   Median : NA  
#>                     Mean   :NaN   Mean   :NaN     Mean   :NaN         Mean   :NaN    Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN  
#>                     3rd Qu.: NA   3rd Qu.: NA     3rd Qu.: NA         3rd Qu.: NA    3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA  
#>                     Max.   : NA   Max.   : NA     Max.   : NA         Max.   : NA    Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA  
#>  chao_evenness   hillebrand      total       appearance  disappearance    jaccard         horn          chao     total_abundance_int    shannon       simpson     inv_simpson 
#>  Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA         Min.   : NA   Min.   : NA   Min.   : NA  
#>  1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA         1st Qu.: NA   1st Qu.: NA   1st Qu.: NA  
#>  Median : NA   Median : NA   Median : NA   Median : NA   Median : NA   Median : NA   Median : NA   Median : NA   Median : NA         Median : NA   Median : NA   Median : NA  
#>  Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN         Mean   :NaN   Mean   :NaN   Mean   :NaN  
#>  3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA         3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA  
#>  Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA         Max.   : NA   Max.   : NA   Max.   : NA  
#>     evenness    jaccard_dis    nestedness     turnover     first_year    last_year        span    
#>  Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA  
#>  1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA  
#>  Median : NA   Median : NA   Median : NA   Median : NA   Median : NA   Median : NA   Median : NA  
#>  Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN  
#>  3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA  
#>  Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA

Tremendous trends in abundance:

high_abundance_slope <- slope_df %>%
  arrange(desc(abs(total_abundance)))
high_abundance_slope
#> # A tibble: 2,678 × 22
#>    siteid total_abundance log_total_abundance species_nb log_species_nb chao_richness chao_shannon chao_simpson chao_evenness  jaccard     horn     chao hillebrand    total appearance
#>    <chr>            <dbl>               <dbl>      <dbl>          <dbl>         <dbl>        <dbl>        <dbl>         <dbl>    <dbl>    <dbl>    <dbl>      <dbl>    <dbl>      <dbl>
#>  1 S11735           -756.              -0.971      0.245          1.31         0.157      -0.0149      -0.0144        0.0142   5.32e-4  1.74e-4 -1.40e-2  -0.00328   0.00451    0.00454
#>  2 S6466            -472.              -7.44       0.309          1.38         0.492       0.199        0.127         0.0529  -3.36e-2 -5.56e-2 -4.77e-2  -0.000754  0.00158    0.00192
#>  3 S5180             398.               3.91       0.207          0.822        0.0876     -0.00311     -0.00540      -0.00641 -5.81e-3 -6.29e-3  5.60e-4  -0.00219   0.00323    0.00351
#>  4 S5812            -368.               1.23       0.393          1.56         0.182       0.0284       0.0163       -0.225   -1.15e-3 -7.21e-3 -3.08e-3  -0.00586   0.00742    0.00739
#>  5 S3863             330.              12.6        0.206          1.53        -0.0549     -0.119       -0.0687       -0.0509  -2.94e-2 -4.84e-3 -1.47e-3  -0.0175    0.0156     0.0138 
#>  6 S1787             277.              29.3        0.337          2.69        -0.0560     -0.0376      -0.0278       -0.00917 -4.65e-3  3.29e-3 -1.17e-2   0.0112   -0.0117    -0.00833
#>  7 S8577             256.               4.78       0.168          0.800       -0.0689     -0.0741      -0.0543       -0.0237  -3.41e-3 -3.93e-4  3.02e-3   0.000423  0.00174    0.00328
#>  8 S4501             255.              29.4        0.582         10.8          0.270       0.0904       0.0269        0.0462  -6.28e-4 -1.06e-2 -2.09e-2   0.00610   0.00559    0.0166 
#>  9 S8158            -249.              -6.37       0.108          0.562        0.517       0.305        0.223         0.0962  -2.94e-2 -2.66e-2  2.00e-5   0.000588 -0.00129   -0.00234
#> 10 S6074             245.               2.37       0.198          0.908        0.0963      0.0740       0.0659        0.0228   6.29e-3 -3.99e-3 -2.55e-4  -0.00519   0.00487    0.00428
#> # … with 2,668 more rows, and 7 more variables: disappearance <dbl>, evenness <dbl>, shannon <dbl>, simpson <dbl>, jaccard_dis <dbl>, nestedness <dbl>, turnover <dbl>
site_desc_slope_total_abundance <-
  high_abundance_slope %>%
  .[["siteid"]]

The most tremendous temporal trends in abundance are located in
Illinois and France

ti <- filtered_dataset$location %>%
  filter(siteid %in% site_desc_slope_total_abundance[1:20]) %>%
  st_as_sf(coords = c("longitude", "latitude"), crs = 2154)

sp::plot(ne_countries())
sp::plot(st_geometry(ti), add = TRUE, col = "red", pch = 19)

filtered_dataset$abun_rich_op %>%
  filter(siteid %in% site_desc_slope_total_abundance[1:20]) %>%
  group_by(siteid) %>%
  summarise(across(c("total_abundance", "species_nb"), list(median = median, mean = mean)))
#> # A tibble: 20 × 5
#>    siteid total_abundance_median total_abundance_mean species_nb_median species_nb_mean
#>  * <chr>                   <dbl>                <dbl>             <dbl>           <dbl>
#>  1 S11735                  5948                17416.              18             18.7 
#>  2 S1507                     86                  698.               2             10.4 
#>  3 S1787                   2790.                2498.              16.5           15.7 
#>  4 S2595                    882                 1306.              15             15.9 
#>  5 S3445                   1791                 2112.               5.5            5.31
#>  6 S3863                   1827                 3353.              13             13.6 
#>  7 S3883                    755                 1842.              10              9.46
#>  8 S4325                    644.                 899.               5.5            5.58
#>  9 S4501                   2111                 1924.              16             14.9 
#> 10 S4755                   2517                 2335.              20             20.1 
#> 11 S5180                   4335                 8431.              23             23.3 
#> 12 S5392                   1416                 4011.              19.5           18.9 
#> 13 S5812                   4968                12216.              24.5           25.2 
#> 14 S5987                   1470                 1630.              18             14.2 
#> 15 S6074                   5908.                9804.              23.5           23.1 
#> 16 S6466                   3166                 6465.              27             25.7 
#> 17 S8158                   3610                 5140.              24.5           24.7 
#> 18 S8577                   4474                 7435.              28             26.4 
#> 19 S8852                   4126                 5016.              25             22.8 
#> 20 S9031                   3078                11660.              22             22.2

p_high_slope_abun <- map(site_desc_slope_total_abundance[1:20],
  ~plot_community_data(
    analysis_dataset %>%
    filter(siteid == .x) %>%
    select(siteid, year, total_abundance),
    y = "total_abundance", x = "year",
  smoothing_method = "lm"))
plot_grid(plotlist = p_high_slope_abun, ncol = 3)

s11735_com <- analysis_dataset %>%
  filter(siteid == site_desc_slope_total_abundance[1]) %>%
  arrange(year) %>%
  select(siteid, year, total_abundance)
s11735_com
#> # A tibble: 25 × 3
#>    siteid  year total_abundance
#>    <chr>  <dbl>           <dbl>
#>  1 S11735  1993              66
#>  2 S11735  1994            1366
#>  3 S11735  1995          109088
#>  4 S11735  1996            9997
#>  5 S11735  1997             512
#>  6 S11735  1998            5948
#>  7 S11735  1999           85163
#>  8 S11735  2000           10785
#>  9 S11735  2001           21720
#> 10 S11735  2002           11719
#> # … with 15 more rows
plot_community_data(s11735_com, y = "total_abundance", x = "year",
  smoothing_method = "lm")

mean(s11735_com$total_abundance) -  sd(s11735_com$total_abundance)
#> [1] -10662.47
mean(s11735_com$total_abundance) +  5 * sd(s11735_com$total_abundance)
#> [1] 157805.5

Petit poissons Notropis atherinoides and Notropis volucellus have several tremendous abundance. Since those bursts happened four times, there have little probability to be error in the data.

s11735_pop <- filtered_dataset$measurement %>%
  filter(siteid == site_desc_slope_total_abundance[1]) %>%
  arrange(year) %>%
  select(siteid, year, species, abundance)

s11735_pop %>%
  group_by(species) %>%
  summarise(abun = sum(abundance)) %>%
  arrange(desc(abun))
#> # A tibble: 52 × 2
#>    species                   abun
#>    <chr>                    <dbl>
#>  1 Notropis atherinoides   418706
#>  2 Notropis volucellus       8400
#>  3 Cyprinella spiloptera     1949
#>  4 Lepomis macrochirus       1659
#>  5 Dorosoma cepedianum        651
#>  6 Aplodinotus grunniens      614
#>  7 Morone chrysops            540
#>  8 Pimephales vigilax         527
#>  9 Cyprinus carpio            427
#> 10 Macrhybopsis aestivalis    386
#> # … with 42 more rows

plot_temporal_population(s11735_pop) +
  theme(legend.position = "bottom")

Second highest temporal trends in abundance:

s6466_com <- analysis_dataset %>%
  filter(siteid == site_desc_slope_total_abundance[2]) %>%
  arrange(year) %>%
  select(siteid, year, total_abundance)
s6466_com
#> # A tibble: 11 × 3
#>    siteid  year total_abundance
#>    <chr>  <dbl>           <dbl>
#>  1 S6466   1981            7844
#>  2 S6466   1983           11938
#>  3 S6466   1987            9521
#>  4 S6466   1988           23997
#>  5 S6466   1989            1558
#>  6 S6466   1990            1674
#>  7 S6466   1991            3166
#>  8 S6466   1992             862
#>  9 S6466   1997            2687
#> 10 S6466   1999            1558
#> 11 S6466   2000            6314
plot_community_data(s6466_com, y = "total_abundance", x = "year",
  smoothing_method = "lm")

s6466_pop <- filtered_dataset$measurement %>%
  filter(siteid == site_desc_slope_total_abundance[2]) %>%
  arrange(year) %>%
  select(siteid, year, species, abundance)

s6466_pop %>%
  group_by(species) %>%
  summarise(abun = sum(abundance)) %>%
  arrange(desc(abun))
#> # A tibble: 55 × 2
#>    species                abun
#>    <chr>                 <dbl>
#>  1 Dorosoma cepedianum   43869
#>  2 Aplodinotus grunniens  8914
#>  3 Notropis wickliffi     8637
#>  4 Alosa chrysochloris    1812
#>  5 Morone chrysops        1538
#>  6 Cyprinus carpio        1272
#>  7 Ictalurus punctatus    1192
#>  8 Notropis atherinoides   745
#>  9 Ictiobus bubalus        559
#> 10 Sander canadensis       465
#> # … with 45 more rows

plot_temporal_population(s6466_pop) +
  theme(legend.position = "bottom")

2 Relationships among temporal trends

2.1 Globally

cor_slope <- slope_df %>%
  select(-siteid) %>%
  cor(., method = "spearman")
cor_slope
#>                     total_abundance log_total_abundance   species_nb log_species_nb chao_richness chao_shannon chao_simpson chao_evenness      jaccard         horn         chao
#> total_abundance         1.000000000         0.830317697  0.277156569    0.250125492    0.01669974 -0.060971568 -0.068607673  -0.046333174 -0.034784436  0.046626092  0.080611195
#> log_total_abundance     0.830317697         1.000000000  0.335207093    0.341531037    0.07832754 -0.004315524 -0.030703462   0.001842386  0.025685809  0.094698554  0.137592452
#> species_nb              0.277156569         0.335207093  1.000000000    0.938150124    0.81502519  0.598119627  0.486160737   0.158532100  0.008257763 -0.012567803 -0.017802771
#> log_species_nb          0.250125492         0.341531037  0.938150124    1.000000000    0.79177426  0.652142713  0.557985838   0.238494063  0.015153040 -0.004767294  0.001695008
#> chao_richness           0.016699737         0.078327544  0.815025195    0.791774260    1.00000000  0.715598168  0.585062363   0.178610995  0.023280895 -0.017989761 -0.051003088
#> chao_shannon           -0.060971568        -0.004315524  0.598119627    0.652142713    0.71559817  1.000000000  0.944296056   0.627245251  0.017092770 -0.041646520 -0.070761250
#> chao_simpson           -0.068607673        -0.030703462  0.486160737    0.557985838    0.58506236  0.944296056  1.000000000   0.673326844  0.001769782 -0.055713026 -0.087567081
#> chao_evenness          -0.046333174         0.001842386  0.158532100    0.238494063    0.17861100  0.627245251  0.673326844   1.000000000  0.069726825  0.032975571  0.039982689
#> jaccard                -0.034784436         0.025685809  0.008257763    0.015153040    0.02328089  0.017092770  0.001769782   0.069726825  1.000000000  0.535041875  0.360518080
#> horn                    0.046626092         0.094698554 -0.012567803   -0.004767294   -0.01798976 -0.041646520 -0.055713026   0.032975571  0.535041875  1.000000000  0.633224505
#> chao                    0.080611195         0.137592452 -0.017802771    0.001695008   -0.05100309 -0.070761250 -0.087567081   0.039982689  0.360518080  0.633224505  1.000000000
#> hillebrand              0.029479915         0.058766270 -0.023380244   -0.022907805   -0.03149248 -0.016465610 -0.027163601   0.063461252  0.272412510  0.335065599  0.629864112
#> total                   0.020626140         0.011421991  0.205181084    0.197786137    0.18733040  0.126450573  0.112381282  -0.042854829 -0.255210721 -0.316495191 -0.589439704
#> appearance              0.161650857         0.192922044  0.687912675    0.691348158    0.58925101  0.453612568  0.384514986   0.090571956 -0.140290491 -0.201907697 -0.347493182
#> disappearance          -0.161078700        -0.222194738 -0.555815437   -0.577843872   -0.45943156 -0.378710413 -0.308701065  -0.143240688 -0.200173498 -0.225399007 -0.445526259
#> evenness               -0.135468004        -0.106903064  0.191546294    0.289979041    0.32190903  0.722795822  0.778140322   0.712515794  0.028733429 -0.023213578 -0.037192382
#> shannon                -0.030652568         0.031284774  0.595797755    0.652755093    0.67989789  0.944628419  0.911228727   0.565547657  0.010823334 -0.044699272 -0.079869118
#> simpson                -0.036762181         0.012861550  0.475730622    0.555069101    0.55478316  0.873409696  0.891031340   0.595960699  0.013305625 -0.042511617 -0.075603823
#> jaccard_dis            -0.028841716        -0.058365304  0.023323286    0.022736545    0.03132624  0.016070172  0.026651172  -0.064111974 -0.272670383 -0.334809565 -0.629675234
#> nestedness              0.004895141        -0.007242935  0.105858491    0.097239453    0.10838694  0.094533541  0.072880587  -0.008434623 -0.075857013 -0.055680929 -0.196458291
#> turnover               -0.022365349        -0.038795778 -0.052263750   -0.048541424   -0.05721969 -0.064561731 -0.031071031  -0.037696468 -0.161988329 -0.230330810 -0.368387266
#>                       hillebrand       total  appearance disappearance     evenness      shannon      simpson  jaccard_dis   nestedness    turnover
#> total_abundance      0.029479915  0.02062614  0.16165086   -0.16107870 -0.135468004 -0.030652568 -0.036762181 -0.028841716  0.004895141 -0.02236535
#> log_total_abundance  0.058766270  0.01142199  0.19292204   -0.22219474 -0.106903064  0.031284774  0.012861550 -0.058365304 -0.007242935 -0.03879578
#> species_nb          -0.023380244  0.20518108  0.68791268   -0.55581544  0.191546294  0.595797755  0.475730622  0.023323286  0.105858491 -0.05226375
#> log_species_nb      -0.022907805  0.19778614  0.69134816   -0.57784387  0.289979041  0.652755093  0.555069101  0.022736545  0.097239453 -0.04854142
#> chao_richness       -0.031492477  0.18733040  0.58925101   -0.45943156  0.321909029  0.679897895  0.554783164  0.031326242  0.108386938 -0.05721969
#> chao_shannon        -0.016465610  0.12645057  0.45361257   -0.37871041  0.722795822  0.944628419  0.873409696  0.016070172  0.094533541 -0.06456173
#> chao_simpson        -0.027163601  0.11238128  0.38451499   -0.30870107  0.778140322  0.911228727  0.891031340  0.026651172  0.072880587 -0.03107103
#> chao_evenness        0.063461252 -0.04285483  0.09057196   -0.14324069  0.712515794  0.565547657  0.595960699 -0.064111974 -0.008434623 -0.03769647
#> jaccard              0.272412510 -0.25521072 -0.14029049   -0.20017350  0.028733429  0.010823334  0.013305625 -0.272670383 -0.075857013 -0.16198833
#> horn                 0.335065599 -0.31649519 -0.20190770   -0.22539901 -0.023213578 -0.044699272 -0.042511617 -0.334809565 -0.055680929 -0.23033081
#> chao                 0.629864112 -0.58943970 -0.34749318   -0.44552626 -0.037192382 -0.079869118 -0.075603823 -0.629675234 -0.196458291 -0.36838727
#> hillebrand           1.000000000 -0.93386273 -0.54773403   -0.64292320  0.047845204 -0.005532024 -0.002397405 -0.999860905 -0.405862238 -0.49276974
#> total               -0.933862733  1.00000000  0.72528811    0.48889404 -0.014337815  0.117869711  0.089656212  0.934006505  0.383373875  0.47223951
#> appearance          -0.547734030  0.72528811  1.00000000   -0.11647853  0.139643932  0.447139132  0.367458109  0.547880499  0.244653628  0.30131921
#> disappearance       -0.642923199  0.48889404 -0.11647853    1.00000000 -0.175733040 -0.387337667 -0.322429697  0.643086266  0.082274006  0.50349434
#> evenness             0.047845204 -0.01433782  0.13964393   -0.17573304  1.000000000  0.784047749  0.866977324 -0.048397766  0.009438247 -0.05364205
#> shannon             -0.005532024  0.11786971  0.44713913   -0.38733767  0.784047749  1.000000000  0.962316417  0.005099809  0.090073607 -0.07159689
#> simpson             -0.002397405  0.08965621  0.36745811   -0.32242970  0.866977324  0.962316417  1.000000000  0.001943911  0.065840971 -0.05210250
#> jaccard_dis         -0.999860905  0.93400650  0.54788050    0.64308627 -0.048397766  0.005099809  0.001943911  1.000000000  0.405983665  0.49279967
#> nestedness          -0.405862238  0.38337388  0.24465363    0.08227401  0.009438247  0.090073607  0.065840971  0.405983665  1.000000000 -0.47739313
#> turnover            -0.492769743  0.47223951  0.30131921    0.50349434 -0.053642052 -0.071596893 -0.052102504  0.492799673 -0.477393127  1.00000000

png(here("doc", "fig", "p_cor_slope.png"))
corrplot::corrplot(cor_slope[c("chao_richness", "species_nb", "log_species_nb",
      "total_abundance", "log_total_abundance"), ], type = "upper", tl.srt = 45, diag = FALSE)
dev.off()
#> X11cairo 
#>        2

png(here("doc", "fig", "p_cor_slope_tot.png"))
corrplot::corrplot(cor_slope, type = "upper", tl.srt = 35, diag = FALSE)
dev.off()
#> X11cairo 
#>        2

corrplot::corrplot(cor_slope[c("chao_richness", "species_nb", "log_species_nb",
      "total_abundance", "log_total_abundance"), ], type = "upper", tl.srt = 45, diag = FALSE)


corrplot::corrplot(cor_slope, type = "upper", tl.srt = 35, diag = FALSE)

ti <- expand.grid(
  resp1 = unique(slope$response),
  resp2 = unique(slope$response)
  ) %>%
  filter(resp2 != resp1) %>%
  filter(
    resp1 %in% c("chao_richness", "species_nb", "log_species_nb",
      "total_abundance", "log_total_abundance")) %>%
  mutate_all(as.character) %>%
  arrange(resp1)

test <- map2(ti$resp1, ti$resp2,
  function(x, y) {
    bi <- slope %>%
      filter(response %in% c(x, y)) %>%
      select(siteid, response, linear_slope) %>%
      pivot_wider(names_from = "response", values_from = "linear_slope")

    return(bi)
  }
)

p_trends_trends <- map(test, function(x) {
  l <- colnames(x)
  x %>%
    ggplot(aes(x = !!sym(l[2]), y = !!sym(l[3]))) +
    geom_point() +
    geom_smooth(method = "loess") +
    labs(
      x = get_var_replacement()[l[2]],
      y = get_var_replacement()[l[3]]
    )
})
names(p_trends_trends) <- map_chr(test, ~colnames(.x)[2])

2.2 Total abundance

plot_grid(
  plotlist = p_trends_trends[names(p_trends_trends) %in% "total_abundance"],
  ncol = 3)
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'

2.3 Species richness

plot_grid(
  plotlist = p_trends_trends[names(p_trends_trends) %in% "species_nb"],
  ncol = 3)
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'

2.4 Log species richness

plot_grid(
  plotlist = p_trends_trends[names(p_trends_trends) %in% "log_species_nb"],
  ncol = 3
)
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'

Trends of total abundance are linked to the one of species richness but not to shannon and simpson, suggesting that the trends of abundance are more link to trends of rare species. Logic: gain and loss of species should concern rare species.
Species richness decrease linked to disappearance, but not increase
Species richness increase linked to appearance, but not decrease

2.5 Turnover vs nestedness:

test <- slope_df %>%
  mutate(
    nes_sup_tu = nestedness > turnover,
    nes_inf_tu = nestedness < turnover,
    rich_inc = species_nb >= 0
  )
sum(test$nes_sup_tu)
#> [1] 1224
sum(test$nes_inf_tu)
#> [1] 1318
test %>%
  tabyl(rich_inc, nes_sup_tu) %>%
  adorn_totals("row") %>%
  adorn_title(placement = "top")
#>           nes_sup_tu     
#>  rich_inc      FALSE TRUE
#>     FALSE        733  499
#>      TRUE        721  725
#>     Total       1454 1224

slope_df_med <- slope_df %>%
  left_join(median_site_rich, by = "siteid")

slope_df_med %>%
  ggplot(aes(
      x = log_species_nb,
      y = nestedness - turnover,
      color = log(med_chao_richness)
      )) +
  geom_point() +
  scale_color_viridis() +
  geom_smooth(method = "loess")
#> `geom_smooth()` using formula 'y ~ x'

median_site_rich <- analysis_dataset %>%
  select(c("siteid", starts_with("chao_"), "total_abundance")) %>%
  group_by(siteid) %>%
  summarise(across(where(is.double), median)) %>%
  rename_with(~paste0("med_", .x), where(is.double))

p_baselga_jaccard <- map(c("nestedness", "turnover"),
  ~slope_df_med %>%
    ggplot(aes(
        x = jaccard_dis,
        y = !!sym(.x),
      color = log(med_chao_richness))) +
    geom_point() +
    scale_color_viridis() +
    geom_abline(intercept = 0, slope = 1) +
    geom_smooth(method = "loess")
)
p_baselga_jaccard[[3]] <-
  slope_df_med %>%
  ggplot(aes(
      x = jaccard_dis,
      y = nestedness - turnover,
      color = log(med_chao_richness))) +
  geom_point() +
  scale_color_viridis() +
  geom_abline(intercept = 0, slope = 1) +
  geom_smooth(method = "loess")

plot_grid(plotlist = p_baselga_jaccard, ncol = 3)
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'

slope_df_med %>%
  ggplot(aes(x = med_chao_richness, y = jaccard_dis)) +
  geom_point()

3 Look at strong turnover trends

Check with Jaccard:

high_jaccard_slope <-
  slope_df %>%
  select(siteid, total_abundance, species_nb, jaccard) %>%
  arrange(desc(jaccard))
high_jaccard_slope
#> # A tibble: 2,678 × 4
#>    siteid total_abundance species_nb jaccard
#>    <chr>            <dbl>      <dbl>   <dbl>
#>  1 S2412           -0.683   0         0.0549
#>  2 S7492            0.205  -4.20e- 2  0.0494
#>  3 S81              5.27   -5.45e- 2  0.0463
#>  4 S1442          -21.1     1.94e- 1  0.0455
#>  5 S1451           -0.175  -2.42e- 2  0.0450
#>  6 S6303           -2.01    0         0.0448
#>  7 S4285            3.35   -2.22e- 2  0.0421
#>  8 S1010           -0.483  -2.38e- 2  0.0420
#>  9 S2613           -6.51   -2.60e- 2  0.0386
#> 10 S869            -3.11   -5.92e-17  0.0376
#> # … with 2,668 more rows
site_desc_jaccard_slope <-
  high_jaccard_slope %>%
  .[["siteid"]]

ti <- filtered_dataset$location %>%
  filter(siteid %in% site_desc_jaccard_slope[1:20]) %>%
  st_as_sf(coords = c("longitude", "latitude"), crs = 2154)

sp::plot(ne_countries())
sp::plot(st_geometry(ti), add = TRUE, col = "red", pch = 19)

filtered_dataset$abun_rich_op %>%
  filter(siteid %in% site_desc_jaccard_slope[1:20]) %>%
  group_by(siteid) %>%
  summarise(across(c("total_abundance", "species_nb"),
      list(median = median, mean = mean))) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "bordered", "hover")) %>%
  scroll_box(width = "100%", height = "400px")

siteid	total_abundance_median	total_abundance_mean	species_nb_median	species_nb_mean
S1010	17.600	23.900000	2.0	2.454546
S1151	1214.000	1278.500000	12.0	11.750000
S1442	178.500	278.400000	8.0	8.000000
S1451	17.000	16.200000	2.0	2.500000
S2102	40.350	49.300000	2.0	1.916667
S2412	12.500	14.200000	1.0	1.000000
S2613	103.500	108.834600	3.0	2.900000
S4285	92.900	120.900000	1.0	1.066667
S6013	96.500	105.700000	8.0	8.400000
S6303	40.050	39.600000	1.0	1.000000
S6550	6.650	14.557143	2.0	2.214286
S6868	57.927	66.733364	2.0	1.636364
S7036	27.600	31.958333	2.0	1.666667
S7329	44.000	44.237500	3.0	3.250000
S7421	10.500	9.835294	3.0	2.764706
S7492	6.450	8.025000	2.5	2.333333
S7805	12.950	24.916667	1.5	1.500000
S81	75.000	73.636364	2.0	2.363636
S8567	103.000	128.266667	4.0	4.400000
S869	38.250	40.225000	1.0	1.000000

p_high_jaccard_slope <- map(site_desc_jaccard_slope[1:20],
  ~plot_community_data(
    analysis_dataset %>%
    filter(siteid == .x) %>%
    select(siteid, year, jaccard),
    y = "jaccard", x = "year",
  smoothing_method = "lm"))
plot_grid(plotlist = p_high_jaccard_slope, ncol = 3)

tar_load(vegdist_turnover_c)

Hum, do not understand why community are more similar at the end of the time-steps

p_pop_high_jaccard <-
  map(site_desc_jaccard_slope[1:20],
    ~plot_temporal_population(
      filtered_dataset$measurement %>%
        filter(siteid == .x)
      ) +
    theme(legend.position = "none")
    )
plot_grid(plotlist = p_pop_high_jaccard, ncol = 3)

Lets see communities that becomes more dissimilar:

p_jaccard_low_jaccard <-
  map(site_desc_jaccard_slope[(length(site_desc_jaccard_slope)-20):length(site_desc_jaccard_slope)],
    ~plot_community_data(
      analysis_dataset %>%
        filter(siteid == .x) %>%
        select(siteid, year, jaccard),
      y = "jaccard", x = "year",
      smoothing_method = "lm"))
plot_grid(plotlist = p_jaccard_low_jaccard, ncol = 3)

p_pop_low_jaccard <-
  map(site_desc_jaccard_slope[
    (length(site_desc_jaccard_slope)-20):length(site_desc_jaccard_slope)],
    ~plot_temporal_population(
      filtered_dataset$measurement %>%
        filter(siteid == .x)
      ) +
    theme(legend.position = "none")
    )
plot_grid(plotlist = p_pop_low_jaccard, ncol = 3)

4 Comparison rigal, lm

rigal_trends_df <-  map2_dfr(
  rigal_trends, names(rigal_trends),
  ~.x %>% mutate(variable = .y)
  ) %>%
  select(variable, siteid, linear_slope) %>%
  pivot_wider(names_from = "variable", values_from = "linear_slope") %>%
  mutate(type = "rigal")

slope_comp <- list()
slope_comp$rigal <- rigal_trends_df[order(rigal_trends_df$siteid), ]
slope_comp$simple_lm <- slope_df[order(slope_df$siteid), ]
stopifnot(all(slope_comp$rigal$siteid == slope_comp$simple_lm$siteid))
slope_comp$diff <- tibble(
  siteid = slope_comp$rigal$siteid)
for (i in var_temporal_trends) {
  slope_comp$diff[[i]] <- slope_comp$rigal[[i]] - slope_comp$simple_lm[[i]]
}

All good !

old_par <- par()
par(mfrow = c(3, 2))
for (i in var_temporal_trends) {
  plot(
    slope_comp$rigal[[i]] ~
      slope_comp$simple_lm[[i]],
    ylab = "Rigal (polynomial degree 2)",
    xlab = "Simple LM"
  )
  abline(0, 1)
  title(i)
}

par(old_par)

# Compare Rigal and simple lm for few sites
plot_rigal_raw_data <- function(
  site = NULL,
  response = NULL,
  dataset = NULL,
  rigal_trends = NULL
  ) {

  raw_data <- dataset %>%
    filter(siteid == site)
  rigal_resp <- rigal_trends[[response]]
  rigal_resp_site <- rigal_resp[rigal_resp$siteid == site, ]
  plot(raw_data[[response]]~raw_data[["year"]], ylab = response, xlab = "Year")
  abline(rigal_resp_site$intercept, rigal_resp_site$linear_slope, col
  = "red")
  abline(lm(raw_data[[response]]~raw_data[["year"]]), col = "green")
  title(main = site)
}

plot_rigal_raw_data(site = "S2672", response = "log_total_abundance", dataset =
analysis_dataset, rigal_trends = rigal_trends)

slope_comp$diff <- slope_comp$diff %>%
  arrange(desc(abs(log_total_abundance)))

old_par <- par()
par(mfrow = c(3, 2))
for (i in slope_comp$diff[1:6,]$siteid) {
  plot_rigal_raw_data(
    site = i,
    response = "log_total_abundance",
    dataset = analysis_dataset,
    rigal_trends = rigal_trends
  )
}

par(old_par)

5 Classification of temporal trends

Classification:

classification_station <- map2(rigal_trends, names(rigal_trends),
  ~.x %>%
  ggplot(aes(x = shape_class)) +
  geom_bar() +
  labs(x = "Trajectory classification", y = "Nb station",
    title = get_var_replacement()[.y]) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
)


plot_grid(plotlist = classification_station)

p_classif_abun_rich <- plot_grid(
  plotlist = classification_station[
    names(classification_station)[str_detect(names(classification_station),
  "abundance|richness|species")]])
save_plot(
  here("doc", "fig", "classif_trends_abun_rich.png"),
  p_classif_abun_rich
)
p_classif_hill_evenness <- plot_grid(
  plotlist = classification_station[
    names(classification_station)[str_detect(names(classification_station),
  "shannon|simpson|evennes")]])
save_plot(
  here("doc", "fig", "classif_trends_hill_evenness.png"),
  p_classif_hill_evenness
)
p_classif_turnover <- plot_grid(
  plotlist = classification_station[
    names(classification_station)[names(classification_station) %in%
      c("jaccard", "horn", "chao", "hillebrand", "total", "appearance",
        "disappearance")]])
save_plot(
  here("doc", "fig", "classif_trends_classif_hill_evenness.png"),
  p_classif_turnover
)

ti <- map_dfr(rigal_trends, ~tabyl_df(x = .x, group = "direction"),
  .id = "response"
)
ti %>%
  filter(direction != "Total") %>%
  select(response, direction, percent) %>%
  pivot_wider(names_from = "direction", values_from = "percent") %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "bordered", "hover"))

response	stable	increase	decrease
total_abundance	75.8%	14.3%	9.9%
log_total_abundance	73.3%	16.5%	10.2%
species_nb	77.7%	13.2%	7.4%
log_species_nb	76.2%	13.2%	6.6%
chao_richness	79.8%	10.6%	8.2%
chao_shannon	76.5%	12.1%	7.7%
chao_simpson	78.4%	10.5%	7.4%
chao_evenness	82.8%	7.5%	5.9%
jaccard	76.0%	3.5%	20.4%
horn	75.9%	4.4%	17.6%
chao	75.7%	4.0%	17.4%
hillebrand	76.8%	2.8%	17.5%
total	73.7%	18.3%	3.0%
appearance	70.2%	15.6%	3.8%
disappearance	64.0%	9.5%	3.0%
evenness	78.8%	9.4%	8.1%
shannon	75.7%	12.2%	8.3%
simpson	77.9%	10.8%	7.7%
jaccard_dis	74.7%	17.4%	2.8%
nestedness	80.3%	9.3%	5.3%
turnover	57.4%	8.3%	1.7%

Linear trends:

#3.3   Log-linear model:logYi=α+βXi+εiIn

#the log-linear model, the literal interpretation of the estimated
#coefficientˆβis that a one-unitincrease inXwill produce an expected increase in
#logYofˆβunits. 
#In terms ofYitself, this
#meansthat  the  expected  value  ofYis  multiplied  byeˆβ.   So  in  terms  of
#effects  of  changes  inXonY(unlogged):•Each 1-unit increase inXmultiplies the
#expected value ofYbyeˆβ.•To compute the effects onYof another change inXthan an
#increase of one unit, call thischangec, we need to includecin the exponent.  The
#effect of ac-unit increase inXis tomultiply the expected value ofYbyecˆβ. So the
#effect for a 5-unit increase inXwould bee5ˆβ.•For small values ofˆβ,
#approximatelyeˆβ≈1+ˆβ. We can use this for the following approxima-tion for a
#quick interpretation of the coefficients:  100·ˆβis the expected percentage
#changeinYfor a unit increase inX.  For instance forˆβ=.06,e.06≈1.06, so a 1-unit
#change inXcorresponds to (approximately) an expected increase inYof 6%

#https://data.library.virginia.edu/interpreting-log-transformations-in-a-linear-model/
hist_linear_plot <- map2(rigal_trends, names(rigal_trends),
  ~.x %>%
  ggplot(aes(x = linear_slope)) +
  geom_histogram() +
  labs(x = "Linear slope", y = "Number of station",
    title = get_var_replacement()[.y]) +
  ylim(c(0, 1000))
)
plot_grid(plotlist = hist_linear_plot)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

type_response <- list(
  abun_rich = str_detect(names(classification_station),
    "abundance|richness|species"),
  hill_evenness = str_detect(names(classification_station),
    "shannon|simpson|evennes"),
  turnover = names(classification_station) %in%
    c("jaccard", "horn", "chao", "hillebrand", "total", "appearance",
      "disappearance")
)

map2(type_response, names(type_response),
  function(x, y) {
  p <- plot_grid(
    plotlist = hist_linear_plot[
      names(hist_linear_plot)[x]
      ]
  )
  save_plot(
    here("doc", "fig", paste0("hist_linear_slope_", y, ".png")),
    p
  )

  }
)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> $abun_rich
#> NULL
#> 
#> $hill_evenness
#> NULL
#> 
#> $turnover
#> NULL

6 Spatialization of temporal trends

rigal_trends_df_loc <-
  rigal_trends_df %>%
  left_join(
    select(filtered_dataset$location, siteid, ecoregion, country),
    by = c("siteid"))

No obvious geographic patterns

rigal_trends_df_loc %>%
  ggplot(aes(x = country, y = log_species_nb, fill = ecoregion)) +
  geom_boxplot()


rigal_trends_df_loc %>%
  ggplot(aes(x = country, y = log_total_abundance, fill = ecoregion)) +
  geom_boxplot()


rigal_trends_df_loc %>%
  ggplot(aes(x = country, y = total, fill = ecoregion)) +
  geom_boxplot()


rigal_trends_df_loc %>%
  ggplot(aes(x = country, y = jaccard, fill = ecoregion)) +
  geom_boxplot()

rigal_trends_df_loc <- filtered_dataset$location %>%
  left_join(rigal_trends_df, by = "siteid") %>%
  st_as_sf(coords = c("longitude", "latitude"),
  crs = 4326)

library(colorspace)
div_pal <- colorspace::sequential_hcl(11, "plasma")
world <- ne_countries(returnclass = "sf")

ggplot(data = world) +
    geom_sf() +
    geom_sf(data = rigal_trends_df_loc,
      aes(color = log_species_nb), shape = 1) +
    scale_color_continuous(type = "viridis") +
    theme(legend.position = "bottom")

rigal_trends_df %>%
  summarise(across(where(is.numeric), median, na.rm = T)) %>%
  unlist()
#>     total_abundance log_total_abundance          species_nb      log_species_nb       chao_richness        chao_shannon        chao_simpson       chao_evenness             jaccard 
#>        0.1129827405        0.0049806811        0.0020269075        0.0004978780        0.0002926352        0.0007389185        0.0007189205        0.0000000000       -0.0045653921 
#>                horn                chao          hillebrand               total          appearance       disappearance            evenness             shannon             simpson 
#>       -0.0022471349       -0.0015502502       -0.0037509723        0.0032257174        0.0009691423        0.0000000000        0.0000000000        0.0004650060        0.0001071122 
#>         jaccard_dis          nestedness            turnover 
#>        0.0037509723        0.0002880739        0.0000000000

median(exp(rigal_trends_df$log_total_abundance) - 1)  * 100
#> [1] 0.4993106
sd(exp(rigal_trends_df$log_total_abundance) - 1) * 100
#> [1] 6.675692
median(exp(na.omit(rigal_trends_df$log_species_nb))  - 1) * 100
#> [1] 0.04980022
sd(exp(na.omit(rigal_trends_df$log_species_nb)) - 1) * 100
#> [1] 2.65351

7 Sensibility analysis of temporal trends

tar_load(c(rigal_trends_avg3y))
names(rigal_trends_avg3y) <- var_temporal_trends

tar_load(c(analysis_dataset_avg3y, measurement_avg3y, measurement))
tar_load(c(turnover_avg3y_c))
tar_load(c(vegdist_turnover_c))
tar_load(c(hillnb_avg3y))

slope_avg3y <- map_dfr(rigal_trends_avg3y[
  ! names(rigal_trends_avg3y) %in% c("jaccard_dis", "nestedness", "turnover")], ~.x %>%
  select(siteid, linear_slope),
.id = "response"
)

Number of sites that are in the different datasets:

u_site <- map(list(classic = slope, avg_3y = slope_avg3y), ~unique(.x$siteid))
sum(u_site[[1]] %in% u_site[[2]])
#> [1] 2678
sum(u_site[[2]] %in% u_site[[1]])
#> [1] 2678

comparison_baseline <- rbind(
  slope %>%
    mutate(type = "classic") %>%
    filter(siteid %in% u_site[["avg_3y"]]),
  slope_avg3y %>%
    mutate(type = "avg_3y")
  ) %>%
pivot_wider(names_from = type, values_from = linear_slope)

It is not too bad, but we can see some differences:
- slope negatively biaised for log_species_nb in 3 years baseline
- Bunch of trends in appearance/disappearance trends disappears with
  3 years baseline
- But note that in three years baseline, there are three points less data

comparison_baseline %>%
  ggplot(aes(x = classic, y = avg_3y)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) +
  facet_wrap(~response, scale = "free")

7.1 Reproducibility

Reproducibility receipt

## datetime
Sys.time()
#> [1] "2022-01-13 17:32:04 CST"

## repository
if(requireNamespace('git2r', quietly = TRUE)) {
  git2r::repository()
} else {
  c(
    system2("git", args = c("log", "--name-status", "-1"), stdout = TRUE),
    system2("git", args = c("remote", "-v"), stdout = TRUE)
  )
}
#> Local:    main /home/alain/Documents/post-these/isu/RivFishTimeBiodiversityFacets
#> Head:     [cd28efe] 2022-01-13: build basin metacom matrix + bray distance

## session info
sessionInfo()
#> R version 4.0.5 (2021-03-31)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Debian GNU/Linux 10 (buster)
#> 
#> Matrix products: default
#> BLAS:   /home/alain/.Renv/versions/4.0.5/lib/R/lib/libRblas.so
#> LAPACK: /home/alain/.Renv/versions/4.0.5/lib/R/lib/libRlapack.so
#> 
#> locale:
#>  [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C               LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8     LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
#>  [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] viridis_0.5.1           viridisLite_0.3.0       colorspace_2.0-0        betapart_1.5.4          untb_1.7-4              future_1.21.0           vegan_2.5-7            
#>  [8] lattice_0.20-41         permute_0.9-5           codyn_2.0.5             janitor_2.1.0           cowplot_1.1.1           rnaturalearthdata_0.1.0 rnaturalearth_0.1.0    
#> [15] sf_1.0-4                scales_1.1.1            kableExtra_1.3.1        here_1.0.1              lubridate_1.7.9.2       magrittr_2.0.1          forcats_0.5.1          
#> [22] stringr_1.4.0           dplyr_1.0.4             purrr_0.3.4             readr_2.1.1             tidyr_1.1.2             tibble_3.1.6            ggplot2_3.3.3          
#> [29] tidyverse_1.3.0         tarchetypes_0.3.2       targets_0.8.1           conflicted_1.1.0        rmarkdown_2.11          nvimcom_0.9-123        
#> 
#> loaded via a namespace (and not attached):
#>   [1] readxl_1.3.1       snow_0.4-4         backports_1.2.1    fastmatch_1.1-3    corrplot_0.84      igraph_1.2.6       sp_1.4-5           splines_4.0.5      gmp_0.6-2.1       
#>  [10] listenv_0.8.0      digest_0.6.27      foreach_1.5.1      htmltools_0.5.1.1  fansi_0.5.0        memoise_2.0.0      cluster_2.1.1      tzdb_0.2.0         globals_0.14.0    
#>  [19] modelr_0.1.8       rvest_0.3.6        haven_2.3.1        xfun_0.28          rgdal_1.5-28       callr_3.7.0        crayon_1.4.2       jsonlite_1.7.2     iterators_1.0.13  
#>  [28] ape_5.5            glue_1.5.1         gtable_0.3.0       webshot_0.5.2      abind_1.4-5        DBI_1.1.1          Rcpp_1.0.6         magic_1.5-9        units_0.6-7       
#>  [37] httr_1.4.2         wk_0.5.0           ellipsis_0.3.2     pkgconfig_2.0.3    partitions_1.10-4  farver_2.0.3       sass_0.3.1         dbplyr_2.1.0       utf8_1.2.2        
#>  [46] tidyselect_1.1.1   labeling_0.4.2     rlang_0.4.12       polynom_1.4-0      munsell_0.5.0      cellranger_1.1.0   tools_4.0.5        cachem_1.0.4       cli_3.1.0         
#>  [55] generics_0.1.0     broom_0.7.4        mathjaxr_1.4-0     evaluate_0.14      geometry_0.4.5     fastmap_1.1.0      yaml_2.2.1         processx_3.5.2     knitr_1.36        
#>  [64] fs_1.5.1           s2_1.0.7           nlme_3.1-152       xml2_1.3.2         compiler_4.0.5     rstudioapi_0.13    e1071_1.7-4        reprex_1.0.0       bslib_0.2.4       
#>  [73] stringi_1.7.6      highr_0.9          ps_1.6.0           Brobdingnag_1.2-6  rgeos_0.5-5        Matrix_1.3-2       classInt_0.4-3     vctrs_0.3.8        pillar_1.6.4      
#>  [82] lifecycle_1.0.1    jquerylib_0.1.3    data.table_1.13.6  R6_2.5.1           bookdown_0.24      gridExtra_2.3      KernSmooth_2.23-18 parallelly_1.23.0  codetools_0.2-18  
#>  [91] rcdd_1.5           MASS_7.3-53.1      assertthat_0.2.1   picante_1.8.2      rprojroot_2.0.2    withr_2.4.3        mgcv_1.8-34        parallel_4.0.5     hms_1.1.1         
#> [100] doSNOW_1.0.19      grid_4.0.5         class_7.3-18       snakecase_0.11.0   git2r_0.29.0       itertools_0.1-3

Temporal trends

Alain Danet

13 janvier, 2022

1 Check trends

2 Relationships among temporal trends

2.1 Globally

2.2 Total abundance

2.3 Species richness

2.4 Log species richness

2.5 Turnover vs nestedness:

3 Look at strong turnover trends

4 Comparison rigal, lm

5 Classification of temporal trends

6 Spatialization of temporal trends

7 Sensibility analysis of temporal trends

7.1 Reproducibility